library(tidyverse) # for data wrangling etc
library(rstanarm) # for fitting models in STAN
library(cmdstanr) # for cmdstan
library(brms) # for fitting models in STAN
library(standist) # for exploring distributions
library(coda) # for diagnostics
library(bayesplot) # for diagnostics
library(ggmcmc) # for MCMC diagnostics
library(DHARMa) # for residual diagnostics
library(rstan) # for interfacing with STAN
library(emmeans) # for marginal means etc
library(broom) # for tidying outputs
library(tidybayes) # for more tidying outputs
library(HDInterval) # for HPD intervals
library(ggeffects) # for partial plots
library(broom.mixed) # for summarising models
library(posterior) # for posterior draws
library(ggeffects) # for partial effects plots
library(patchwork) # for multi-panel figures
library(bayestestR) # for ROPE
library(see) # for some plots
library(easystats) # framework for stats, modelling and visualisation
library(geoR) # framework for stats, modelling and visualisation
library(modelsummary) # for data and model summaries
theme_set(theme_grey()) # put the default ggplot theme back
source("helperFunctions.R")Bayesian GLM Part7
1 Preparations
Load the necessary libraries
2 Scenario
Mitchell et al. (2019) presented a data set in which male guppy fish where placed individually into long, narrow and shallow tanks that predominately only permitted the fish to swim in two dimensions. The activity rate (cumulative distance swam) of each guppy fish was remotely recorded via tracking software during 20 minute trials conducted each day for 14 days. The order in which individual guppy fish were recorded each day was randomised and the time of day of each trial was recorded.
The researchers were interested in whether the the guppy fish become habituated to the experimental conditions and altered their activity patterns over time.
The data are in the file mitchell.csv in the data folder.
| ID | TOD | Day | Dist_moved |
|---|---|---|---|
| A39 | 0.39 | 7 | 3439.93 |
| A3 | 0.39 | 7 | 2795.08 |
| A13 | 0.39 | 7 | 5258.79 |
| ... | ... |
| ID: | Individual guppy fish ID - Random variable |
| TOD: | Time of day of observation as a proportion - Predictor variable |
| Day: | Day post entry into tank - Predictor variable |
| Dist_moved: | Total distance (cm) moved in 20 min trial - Response variable |
3 Read in the data
Rows: 1626 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (3): ID, Date, StartDate
dbl (5): Mass, Day, TOD, Dist_moved, Dist
time (1): TimeOfDay
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Rows: 1,626
Columns: 9
$ ID <chr> "A39", "A3", "A13", "A23", "A59", "A69", "A41", "A9", "A35"…
$ TimeOfDay <time> 09:28:32, 09:28:32, 09:28:32, 09:28:32, 09:28:32, 09:28:32…
$ Date <chr> "6/03/2018", "6/03/2018", "6/03/2018", "6/03/2018", "6/03/2…
$ StartDate <chr> "27-Feb-18", "27-Feb-18", "27-Feb-18", "27-Feb-18", "27-Feb…
$ Mass <dbl> 0.066, 0.072, 0.105, 0.080, 0.077, 0.092, 0.073, 0.075, 0.0…
$ Day <dbl> 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,…
$ TOD <dbl> 0.39, 0.39, 0.39, 0.39, 0.39, 0.39, 0.39, 0.42, 0.42, 0.42,…
$ Dist_moved <dbl> 3439.930, 2795.080, 5258.790, 3402.740, 7672.660, 2400.660,…
$ Dist <dbl> 58.65092, 52.86852, 72.51752, 58.33301, 87.59372, 48.99653,…
# A tibble: 6 × 9
ID TimeOfDay Date StartDate Mass Day TOD Dist_moved Dist
<chr> <time> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 A39 09:28:32 6/03/2018 27-Feb-18 0.066 7 0.39 3440. 58.7
2 A3 09:28:32 6/03/2018 27-Feb-18 0.072 7 0.39 2795. 52.9
3 A13 09:28:32 6/03/2018 27-Feb-18 0.105 7 0.39 5259. 72.5
4 A23 09:28:32 6/03/2018 27-Feb-18 0.08 7 0.39 3403. 58.3
5 A59 09:28:32 6/03/2018 27-Feb-18 0.077 7 0.39 7673. 87.6
6 A69 09:28:32 6/03/2018 27-Feb-18 0.092 7 0.39 2401. 49.0
spc_tbl_ [1,626 × 9] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
$ ID : chr [1:1626] "A39" "A3" "A13" "A23" ...
$ TimeOfDay : 'hms' num [1:1626] 09:28:32 09:28:32 09:28:32 09:28:32 ...
..- attr(*, "units")= chr "secs"
$ Date : chr [1:1626] "6/03/2018" "6/03/2018" "6/03/2018" "6/03/2018" ...
$ StartDate : chr [1:1626] "27-Feb-18" "27-Feb-18" "27-Feb-18" "27-Feb-18" ...
$ Mass : num [1:1626] 0.066 0.072 0.105 0.08 0.077 0.092 0.073 0.075 0.09 0.077 ...
$ Day : num [1:1626] 7 7 7 7 7 7 7 7 7 7 ...
$ TOD : num [1:1626] 0.39 0.39 0.39 0.39 0.39 0.39 0.39 0.42 0.42 0.42 ...
$ Dist_moved: num [1:1626] 3440 2795 5259 3403 7673 ...
$ Dist : num [1:1626] 58.7 52.9 72.5 58.3 87.6 ...
- attr(*, "spec")=
.. cols(
.. ID = col_character(),
.. TimeOfDay = col_time(format = ""),
.. Date = col_character(),
.. StartDate = col_character(),
.. Mass = col_double(),
.. Day = col_double(),
.. TOD = col_double(),
.. Dist_moved = col_double(),
.. Dist = col_double()
.. )
- attr(*, "problems")=<externalptr>
mitchell (1626 rows and 9 variables, 9 shown)
ID | Name | Type | Missings | Values | N
---+------------+-----------+----------+-------------------+------------
1 | ID | character | 0 (0.0%) | A1 | 14 ( 0.9%)
| | | | A11 | 14 ( 0.9%)
| | | | A13 | 14 ( 0.9%)
| | | | A15 | 14 ( 0.9%)
| | | | A17 | 14 ( 0.9%)
| | | | A19 | 14 ( 0.9%)
| | | | A21 | 14 ( 0.9%)
| | | | A23 | 14 ( 0.9%)
| | | | A25 | 14 ( 0.9%)
| | | | A27 | 14 ( 0.9%)
| | | | (...) |
---+------------+-----------+----------+-------------------+------------
2 | TimeOfDay | numeric | 0 (0.0%) | 07:57:31 | 7 ( 0.4%)
| | | | 08:08:00 | 7 ( 0.4%)
| | | | 08:16:18 | 2 ( 0.1%)
| | | | 08:17:16 | 8 ( 0.5%)
| | | | 08:17:56 | 5 ( 0.3%)
| | | | 08:18:36 | 7 ( 0.4%)
| | | | 08:20:22 | 6 ( 0.4%)
| | | | 08:23:28 | 8 ( 0.5%)
| | | | 08:26:07 | 5 ( 0.3%)
| | | | 08:30:43 | 10 ( 0.6%)
| | | | (...) |
---+------------+-----------+----------+-------------------+------------
3 | Date | character | 0 (0.0%) | 1/04/2018 | 40 ( 2.5%)
| | | | 1/05/2018 | 40 ( 2.5%)
| | | | 10/03/2018 | 39 ( 2.4%)
| | | | 10/04/2018 | 40 ( 2.5%)
| | | | 11/03/2018 | 39 ( 2.4%)
| | | | 11/04/2018 | 40 ( 2.5%)
| | | | 12/03/2018 | 39 ( 2.4%)
| | | | 12/04/2018 | 40 ( 2.5%)
| | | | 13/03/2018 | 39 ( 2.4%)
| | | | 14/03/2018 | 39 ( 2.4%)
| | | | (...) |
---+------------+-----------+----------+-------------------+------------
4 | StartDate | character | 0 (0.0%) | 17-Apr-18 | 560 (34.4%)
| | | | 23-Mar-18 | 520 (32.0%)
| | | | 27-Feb-18 | 546 (33.6%)
---+------------+-----------+----------+-------------------+------------
5 | Mass | numeric | 0 (0.0%) | [0.05, 0.12] | 1626
---+------------+-----------+----------+-------------------+------------
6 | Day | numeric | 0 (0.0%) | [7, 20] | 1626
---+------------+-----------+----------+-------------------+------------
7 | TOD | numeric | 0 (0.0%) | [0.37, 0.56] | 1626
---+------------+-----------+----------+-------------------+------------
8 | Dist_moved | numeric | 0 (0.0%) | [101.13, 11833.9] | 1626
---+------------+-----------+----------+-------------------+------------
9 | Dist | numeric | 0 (0.0%) | [10.06, 108.78] | 1626
------------------------------------------------------------------------
Warning: These variables were omitted because they include more than 50 levels:
ID.
| Unique | Missing Pct. | Mean | SD | Min | Median | Max | Histogram | |
|---|---|---|---|---|---|---|---|---|
| Mass | 51 | 0 | 0.1 | 0.0 | 0.0 | 0.1 | 0.1 | |
| Day | 14 | 0 | 13.5 | 4.1 | 7.0 | 13.0 | 20.0 | |
| TOD | 97 | 0 | 0.5 | 0.0 | 0.4 | 0.4 | 0.6 | |
| Dist_moved | 1625 | 0 | 3917.0 | 2141.8 | 101.1 | 3811.0 | 11833.9 | |
| Dist | 1625 | 0 | 59.9 | 18.1 | 10.1 | 61.7 | 108.8 | |
| N | % | |||||||
| Date | 1/04/2018 | 40 | 2.5 | |||||
| 1/05/2018 | 40 | 2.5 | ||||||
| 10/03/2018 | 39 | 2.4 | ||||||
| 10/04/2018 | 40 | 2.5 | ||||||
| 11/03/2018 | 39 | 2.4 | ||||||
| 11/04/2018 | 40 | 2.5 | ||||||
| 12/03/2018 | 39 | 2.4 | ||||||
| 12/04/2018 | 40 | 2.5 | ||||||
| 13/03/2018 | 39 | 2.4 | ||||||
| 14/03/2018 | 39 | 2.4 | ||||||
| 15/03/2018 | 39 | 2.4 | ||||||
| 16/03/2018 | 39 | 2.4 | ||||||
| 17/03/2018 | 39 | 2.4 | ||||||
| 18/03/2018 | 39 | 2.4 | ||||||
| 19/03/2018 | 39 | 2.4 | ||||||
| 2/04/2018 | 40 | 2.5 | ||||||
| 2/05/2018 | 40 | 2.5 | ||||||
| 24/04/2018 | 40 | 2.5 | ||||||
| 25/04/2018 | 40 | 2.5 | ||||||
| 26/04/2018 | 40 | 2.5 | ||||||
| 27/04/2018 | 40 | 2.5 | ||||||
| 28/04/2018 | 40 | 2.5 | ||||||
| 29/04/2018 | 40 | 2.5 | ||||||
| 3/04/2018 | 40 | 2.5 | ||||||
| 3/05/2018 | 40 | 2.5 | ||||||
| 30/03/2018 | 40 | 2.5 | ||||||
| 30/04/2018 | 40 | 2.5 | ||||||
| 31/03/2018 | 40 | 2.5 | ||||||
| 4/04/2018 | 40 | 2.5 | ||||||
| 4/05/2018 | 40 | 2.5 | ||||||
| 5/04/2018 | 40 | 2.5 | ||||||
| 5/05/2018 | 40 | 2.5 | ||||||
| 6/03/2018 | 39 | 2.4 | ||||||
| 6/04/2018 | 40 | 2.5 | ||||||
| 6/05/2018 | 40 | 2.5 | ||||||
| 7/03/2018 | 39 | 2.4 | ||||||
| 7/05/2018 | 40 | 2.5 | ||||||
| 8/03/2018 | 39 | 2.4 | ||||||
| 8/04/2018 | 40 | 2.5 | ||||||
| 9/03/2018 | 39 | 2.4 | ||||||
| 9/04/2018 | 40 | 2.5 | ||||||
| StartDate | 17-Apr-18 | 560 | 34.4 | |||||
| 23-Mar-18 | 520 | 32.0 | ||||||
| 27-Feb-18 | 546 | 33.6 |
Warning: These variables were omitted because they include more than 50 levels:
ID.
| Day | Unique | Missing Pct. | Mean | SD | Min | Median | Max | Histogram | |
|---|---|---|---|---|---|---|---|---|---|
| Mass | 7.0 | 51 | 0 | 0.1 | 0.0 | 0.0 | 0.1 | 0.1 | |
| 8.0 | 51 | 0 | 0.1 | 0.0 | 0.0 | 0.1 | 0.1 | ||
| 9.0 | 51 | 0 | 0.1 | 0.0 | 0.0 | 0.1 | 0.1 | ||
| 10.0 | 51 | 0 | 0.1 | 0.0 | 0.0 | 0.1 | 0.1 | ||
| 11.0 | 51 | 0 | 0.1 | 0.0 | 0.0 | 0.1 | 0.1 | ||
| 12.0 | 51 | 0 | 0.1 | 0.0 | 0.0 | 0.1 | 0.1 | ||
| 13.0 | 51 | 0 | 0.1 | 0.0 | 0.0 | 0.1 | 0.1 | ||
| 14.0 | 51 | 0 | 0.1 | 0.0 | 0.0 | 0.1 | 0.1 | ||
| 15.0 | 45 | 0 | 0.1 | 0.0 | 0.0 | 0.1 | 0.1 | ||
| 16.0 | 51 | 0 | 0.1 | 0.0 | 0.0 | 0.1 | 0.1 | ||
| 17.0 | 51 | 0 | 0.1 | 0.0 | 0.0 | 0.1 | 0.1 | ||
| 18.0 | 51 | 0 | 0.1 | 0.0 | 0.0 | 0.1 | 0.1 | ||
| 19.0 | 51 | 0 | 0.1 | 0.0 | 0.0 | 0.1 | 0.1 | ||
| 20.0 | 51 | 0 | 0.1 | 0.0 | 0.0 | 0.1 | 0.1 | ||
| TOD | 7.0 | 12 | 0 | 0.5 | 0.0 | 0.4 | 0.5 | 0.5 | |
| 8.0 | 11 | 0 | 0.5 | 0.1 | 0.4 | 0.4 | 0.6 | ||
| 9.0 | 18 | 0 | 0.5 | 0.1 | 0.4 | 0.4 | 0.6 | ||
| 10.0 | 18 | 0 | 0.5 | 0.0 | 0.4 | 0.5 | 0.5 | ||
| 11.0 | 19 | 0 | 0.4 | 0.0 | 0.4 | 0.4 | 0.5 | ||
| 12.0 | 19 | 0 | 0.4 | 0.0 | 0.4 | 0.4 | 0.5 | ||
| 13.0 | 20 | 0 | 0.5 | 0.1 | 0.4 | 0.4 | 0.5 | ||
| 14.0 | 21 | 0 | 0.5 | 0.1 | 0.4 | 0.4 | 0.5 | ||
| 15.0 | 13 | 0 | 0.5 | 0.1 | 0.4 | 0.5 | 0.6 | ||
| 16.0 | 17 | 0 | 0.5 | 0.0 | 0.4 | 0.4 | 0.5 | ||
| 17.0 | 17 | 0 | 0.5 | 0.0 | 0.4 | 0.5 | 0.6 | ||
| 18.0 | 20 | 0 | 0.5 | 0.1 | 0.4 | 0.5 | 0.5 | ||
| 19.0 | 17 | 0 | 0.5 | 0.0 | 0.4 | 0.4 | 0.6 | ||
| 20.0 | 18 | 0 | 0.5 | 0.0 | 0.4 | 0.5 | 0.5 | ||
| Dist_moved | 7.0 | 119 | 0 | 3766.9 | 2207.2 | 101.1 | 3402.7 | 9037.0 | |
| 8.0 | 119 | 0 | 3609.6 | 2066.8 | 269.8 | 3589.3 | 8343.5 | ||
| 9.0 | 119 | 0 | 3468.2 | 1996.3 | 251.3 | 3396.0 | 8581.9 | ||
| 10.0 | 119 | 0 | 3525.7 | 1908.6 | 472.4 | 3256.8 | 8343.1 | ||
| 11.0 | 119 | 0 | 3603.1 | 2032.4 | 341.6 | 3389.6 | 8508.7 | ||
| 12.0 | 119 | 0 | 3682.0 | 2081.9 | 577.9 | 3425.5 | 9167.5 | ||
| 13.0 | 119 | 0 | 3918.3 | 2181.0 | 415.2 | 3691.7 | 8840.6 | ||
| 14.0 | 119 | 0 | 4018.1 | 2107.9 | 566.2 | 3928.4 | 8881.7 | ||
| 15.0 | 79 | 0 | 4031.3 | 2212.8 | 685.7 | 3851.9 | 9619.0 | ||
| 16.0 | 119 | 0 | 4199.5 | 2241.9 | 438.1 | 4215.6 | 10238.8 | ||
| 17.0 | 119 | 0 | 4234.3 | 2222.7 | 549.3 | 4107.0 | 9439.1 | ||
| 18.0 | 119 | 0 | 4315.3 | 2128.9 | 563.4 | 4684.7 | 9420.1 | ||
| 19.0 | 119 | 0 | 4238.4 | 2202.5 | 656.7 | 4053.9 | 9186.6 | ||
| 20.0 | 119 | 0 | 4266.2 | 2220.8 | 433.4 | 4454.1 | 11833.9 | ||
| Dist | 7.0 | 119 | 0 | 58.3 | 19.4 | 10.1 | 58.3 | 95.1 | |
| 8.0 | 119 | 0 | 57.3 | 18.1 | 16.4 | 59.9 | 91.3 | ||
| 9.0 | 119 | 0 | 56.2 | 17.7 | 15.9 | 58.3 | 92.6 | ||
| 10.0 | 119 | 0 | 56.9 | 16.9 | 21.7 | 57.1 | 91.3 | ||
| 11.0 | 119 | 0 | 57.3 | 18.0 | 18.5 | 58.2 | 92.2 | ||
| 12.0 | 119 | 0 | 58.1 | 17.6 | 24.0 | 58.5 | 95.7 | ||
| 13.0 | 119 | 0 | 59.8 | 18.6 | 20.4 | 60.8 | 94.0 | ||
| 14.0 | 119 | 0 | 60.9 | 17.7 | 23.8 | 62.7 | 94.2 | ||
| 15.0 | 79 | 0 | 60.8 | 18.4 | 26.2 | 62.1 | 98.1 | ||
| 16.0 | 119 | 0 | 62.2 | 18.4 | 20.9 | 64.9 | 101.2 | ||
| 17.0 | 119 | 0 | 62.6 | 18.0 | 23.4 | 64.1 | 97.2 | ||
| 18.0 | 119 | 0 | 63.3 | 17.6 | 23.7 | 68.4 | 97.1 | ||
| 19.0 | 119 | 0 | 62.6 | 18.0 | 25.6 | 63.7 | 95.8 | ||
| 20.0 | 119 | 0 | 62.8 | 18.1 | 20.8 | 66.7 | 108.8 | ||
| N | % | ||||||||
| Date | 1/04/2018 | 40 | 2.5 | ||||||
| 1/05/2018 | 40 | 2.5 | |||||||
| 10/03/2018 | 39 | 2.4 | |||||||
| 10/04/2018 | 40 | 2.5 | |||||||
| 11/03/2018 | 39 | 2.4 | |||||||
| 11/04/2018 | 40 | 2.5 | |||||||
| 12/03/2018 | 39 | 2.4 | |||||||
| 12/04/2018 | 40 | 2.5 | |||||||
| 13/03/2018 | 39 | 2.4 | |||||||
| 14/03/2018 | 39 | 2.4 | |||||||
| 15/03/2018 | 39 | 2.4 | |||||||
| 16/03/2018 | 39 | 2.4 | |||||||
| 17/03/2018 | 39 | 2.4 | |||||||
| 18/03/2018 | 39 | 2.4 | |||||||
| 19/03/2018 | 39 | 2.4 | |||||||
| 2/04/2018 | 40 | 2.5 | |||||||
| 2/05/2018 | 40 | 2.5 | |||||||
| 24/04/2018 | 40 | 2.5 | |||||||
| 25/04/2018 | 40 | 2.5 | |||||||
| 26/04/2018 | 40 | 2.5 | |||||||
| 27/04/2018 | 40 | 2.5 | |||||||
| 28/04/2018 | 40 | 2.5 | |||||||
| 29/04/2018 | 40 | 2.5 | |||||||
| 3/04/2018 | 40 | 2.5 | |||||||
| 3/05/2018 | 40 | 2.5 | |||||||
| 30/03/2018 | 40 | 2.5 | |||||||
| 30/04/2018 | 40 | 2.5 | |||||||
| 31/03/2018 | 40 | 2.5 | |||||||
| 4/04/2018 | 40 | 2.5 | |||||||
| 4/05/2018 | 40 | 2.5 | |||||||
| 5/04/2018 | 40 | 2.5 | |||||||
| 5/05/2018 | 40 | 2.5 | |||||||
| 6/03/2018 | 39 | 2.4 | |||||||
| 6/04/2018 | 40 | 2.5 | |||||||
| 6/05/2018 | 40 | 2.5 | |||||||
| 7/03/2018 | 39 | 2.4 | |||||||
| 7/05/2018 | 40 | 2.5 | |||||||
| 8/03/2018 | 39 | 2.4 | |||||||
| 8/04/2018 | 40 | 2.5 | |||||||
| 9/03/2018 | 39 | 2.4 | |||||||
| 9/04/2018 | 40 | 2.5 | |||||||
| StartDate | 17-Apr-18 | 560 | 34.4 | ||||||
| 23-Mar-18 | 520 | 32.0 | |||||||
| 27-Feb-18 | 546 | 33.6 |
4 Data preparation
Let start by declaring the categorical variables and random effect as factors.
5 Exploratory data analysis
As the response represents a strictly positive real number (total distance moved), it would make sense to start by considering a Gamma model.
Note, also that the response variable are relatively large distances (as the units of distance used were centimeters). It might be worth dividing these values by 100 or even 1000 during the modelling process.
Model formula: \[ y_i/1000 \sim{} \mathcal{Gamma}(\lambda_i)\\ ln(\lambda_i) =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]
where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of day and time of day on the total distance moved during the trial. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual guppies.
mitchell |> ggplot(aes(y = Dist, x = Day)) +
geom_point() +
scale_y_log10() +
scale_x_log10() +
geom_smooth(method = "lm") +
facet_wrap(~ID)`geom_smooth()` using formula = 'y ~ x'
mitchell |> ggplot(aes(y = Dist, x = TOD)) +
geom_point() +
scale_y_log10() +
scale_x_log10() +
geom_smooth(method = "lm") +
facet_wrap(~ID)`geom_smooth()` using formula = 'y ~ x'
6 Fit the model
In brms, the default priors are designed to be weakly informative. They are chosen to provide moderate regularisation (to help prevent over fitting) and help stabilise the computations.
Unlike rstanarm, brms models must be compiled before they start sampling. For most models, the compilation of the stan code takes around 45 seconds.
mitchell.form <- bf(I(Dist_moved/1000) ~ scale(Day) + scale(TOD),
family=Gamma(link='log'))
options(width=150)
mitchell.form |> get_prior(data = mitchell_sub) prior class coef group resp dpar nlpar lb ub source
(flat) b default
(flat) b scaleDay (vectorized)
(flat) b scaleTOD (vectorized)
student_t(3, 1.3, 2.5) Intercept default
gamma(0.01, 0.01) shape 0 default
The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]
When defining our own priors, we typically do not want them to be scaled.
If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucked out of thin air):
- \(\beta_0\): normal centred at 134 with a standard deviation of 65
- mean of 134: since
median(log(mitchell$Dist_moved/1000)) - sd pf 0.7: since
mad(log(mitchell$Dist_moved/1000))
- mean of 134: since
- \(\beta_1\): normal centred at 0 with a standard deviation of 2
- sd of 1: since
mad(log(Dist_moved/1000))/mad(scale(Day))
- sd of 1: since
- \(\sigma\): half-t centred at 0 with a standard deviation of 65 OR
- sd pf 65: since
mad(log(mitchell$Dist_moved/1000))
- sd pf 65: since
- \(\sigma\): gamma with shape parameters of 2 and 1
- \(\beta_0\): normal centred at 1.5 with a standard deviation of 1.5
- mean of 1.5: since
mean(log(mitchell$Dist_moved/1000))ormean(asinh(mitchell$Dist_moved/(2*1))/log(exp(1))) - sd of 1.5: since
sd(log(mitchell$Dist_moved/1000))orsd(asinh(mitchell$Dist_moved/(2*1))/log(exp(1)))
- mean of 1.5: since
mitchell_sub |>
mutate(Dist = Dist) |>
summarise(
Int_mu = median(Dist),
Int_sd = mad(Dist),
b1_sd = mad(Dist) / mad(scale(Day)),
b2_sd = mad(Dist) / mad(scale(TOD))
)# A tibble: 1 × 4
Int_mu Int_sd b1_sd b2_sd
<dbl> <dbl> <dbl> <dbl>
1 61.5 13.4 10.8 9.97
## mitchell |>
## mutate(lDist = log(Dist_moved/1)) |>
## summarise(
## Int_mu = median(lDist),
## Int_sd = mad(lDist),
## b1_sd = mad(lDist)/mad(scale(Day)),
## b2_sd = mad(lDist)/mad(scale(TOD))
## )
## standist::visualize("normal(1.4, 1)", xlim=c(0,20))
## standist::visualize("student_t(3, 0, 1)",
## xlim=c(-10,25))## priors <- prior(normal(1.4, 1), class = 'Intercept') +
## prior(normal(0, 1), class = 'b') +
## prior(student_t(3, 0, 1), class = 'sd')
## prior(gamma(2, 1), class = 'shape')
priors <- prior(normal(60, 25), class = "Intercept") +
prior(normal(0, 15), class = "b") +
## prior(student_t(3, 0, 25), class = 'sd') +
prior(student_t(3, 0, 25), class = "sigma")
## prior(gamma(2, 1), class = 'shape')
## mitchell.form <- bf(Dist_moved ~ scale(Day) + scale(TOD) + (scale(Day) + scale(TOD)|ID),
## ## family=Gamma(link='log'))
## family=gaussian())
mitchell.form <- bf(Dist ~ scale(Day) + scale(TOD),
## family=Gamma(link='log'))
family = gaussian()
)
mitchell.brm2 <- brm(mitchell.form,
data = mitchell_sub,
prior = priors,
sample_prior = "only",
iter = 5000,
warmup = 2500,
chains = 3,
cores = 3,
thin = 10,
refresh = 100,
seed = 123,
control = list(adapt_delta = 0.99, max_treedepth = 20),
backend = "rstan"
)Compiling Stan program...
Start sampling
mitchell.brm2 |>
conditional_effects("Day") |>
plot(points = TRUE) |>
_[[1]] +
scale_y_continuous(trans = scales::pseudo_log_trans())mitchell.brm2 |>
conditional_effects("TOD") |>
plot(points = TRUE) |>
_[[1]] +
scale_y_continuous(trans = scales::pseudo_log_trans())The above seem sufficiently wide whilst at the same time not providing any encouragement for the sampler to wander off into very unsupported areas.
The desired updates require recompiling the model
Compiling Stan program...
Start sampling
7 MCMC sampling diagnostics
The brms package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.
Of these, we will focus on:
- stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).
The chains appear well mixed and very similar
- stan_acf (auto-correlation function): plots the auto-correlation between successive MCMC sample lags for each parameter and each chain
There is no evidence of auto-correlation in the MCMC samples
- stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.
All Rhat values are below 1.05, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
Ratios all very high.
8 Model validation
Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.
bayesplot PPC module:
ppc_bars
ppc_bars_grouped
ppc_boxplot
ppc_dens
ppc_dens_overlay
ppc_dens_overlay_grouped
ppc_ecdf_overlay
ppc_ecdf_overlay_grouped
ppc_error_binned
ppc_error_hist
ppc_error_hist_grouped
ppc_error_scatter
ppc_error_scatter_avg
ppc_error_scatter_avg_grouped
ppc_error_scatter_avg_vs_x
ppc_freqpoly
ppc_freqpoly_grouped
ppc_hist
ppc_intervals
ppc_intervals_grouped
ppc_km_overlay
ppc_km_overlay_grouped
ppc_loo_intervals
ppc_loo_pit_ecdf
ppc_loo_pit_overlay
ppc_loo_pit_qq
ppc_loo_ribbon
ppc_pit_ecdf
ppc_pit_ecdf_grouped
ppc_ribbon
ppc_ribbon_grouped
ppc_rootogram
ppc_scatter
ppc_scatter_avg
ppc_scatter_avg_grouped
ppc_stat
ppc_stat_2d
ppc_stat_freqpoly
ppc_stat_freqpoly_grouped
ppc_stat_grouped
ppc_violin_grouped
- dens_overlay: plots the density distribution of the observed data (black line) overlayed on top of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.
The model draws appear to represent the shape of the observed data reasonably well
- error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot. There is some pattern remaining in these residuals.
Using all posterior draws for ppc type 'error_scatter_avg' by default.
The predictive error seems to be related to the predictor - the model performs poorest at higher recruitments.
- intervals: plots the observed data overlayed on top of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.
Using all posterior draws for ppc type 'intervals' by default.
Whilst the modelled predictions do a reasonable job of representing the observed data, the observed data do appear to be more varied than the model is representing.
The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.
DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components yourself, we can still obtain the simulated residuals from the fitted stan model.
We need to supply:
- simulated (predicted) responses associated with each observation.
- observed values
- fitted (predicted) responses (averaged) associated with each observation
Warning: Argument 'nsamples' is deprecated. Please use argument 'ndraws'
instead.
## mitchell.resids <- createDHARMa(simulatedResponse = t(preds),
## observedResponse = mitchell$Dist_moved/1000,
## fittedPredictedResponse = apply(preds, 2, median),
## integerResponse = TRUE)
mitchell.resids <- createDHARMa(
simulatedResponse = t(preds),
observedResponse = mitchell_sub$Dist,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE
)
mitchell.resids |> plot()
DHARMa nonparametric dispersion test via sd of residuals fitted vs.
simulated
data: simulationOutput
dispersion = 0.49418, p-value = 0.416
alternative hypothesis: two.sided
mitchell.resids <- make_brms_dharma_res(mitchell.brm3, integerResponse = FALSE)
wrap_elements(~ testUniformity(mitchell.resids)) +
wrap_elements(~ plotResiduals(mitchell.resids, form = factor(rep(1, nrow(mitchell_sub))))) +
wrap_elements(~ plotResiduals(mitchell.resids, quantreg = TRUE)) +
wrap_elements(~ testDispersion(mitchell.resids))
Durbin-Watson test
data: simulationOutput$scaledResiduals ~ 1
DW = 2.0389, p-value = 0.9405
alternative hypothesis: true autocorrelation is not 0
## mitchell.resid1 <- mitchell.resids |>
## recalculateResiduals(group = mitchell$Day, aggregateBy = mean)
## ## recalculateResiduals(group = interaction(mitchell$Day, mitchell$ID), aggregateBy = mean)
## mitchell.resid1 |> testTemporalAutocorrelation(time=unique(mitchell$Day))
autocor_check(mitchell_sub, mitchell.brm3, variable = "Day", n.sim = 250)variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
Using as id variables
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
`geom_smooth()` using formula = 'y ~ x'
Conclusions:
- the simulated residuals do suggest that there might be a dispersion issue
- it might be worth exploring either zero-inflation, a negative binomial model, or including a observation-level random effect.
9 Refit model
## prior(gamma(2, 1), class = 'shape')
priors <- prior(normal(60, 10), class = "Intercept") +
prior(normal(0, 10), class = "b") +
## prior(student_t(3, 0, 10), class = 'sd') +
prior(student_t(3, 0, 10), class = "sigma") #+
## prior(lkj_corr_cholesky(1), class = 'cor')
## prior(gamma(0.01, 0.01), class = 'shape')
mitchell.form <- bf(
Dist ~ scale(Day) + scale(TOD) +
ar(time = Day, gr = ID, p = 1, cov = TRUE),
family = gaussian()
)
## family=lognormal())
## family=Gamma(link='log'))
get_prior(mitchell.form, data = mitchell_sub) prior class coef group resp dpar nlpar lb ub
(flat) ar -1 1
(flat) b
(flat) b scaleDay
(flat) b scaleTOD
student_t(3, 61.5, 13.4) Intercept
student_t(3, 0, 13.4) sigma 0
source
default
default
(vectorized)
(vectorized)
default
default
mitchell.brm4 <- brm(mitchell.form,
data = mitchell_sub,
prior = priors,
sample_prior = "yes",
iter = 5000,
warmup = 2500,
chains = 3,
cores = 3,
thin = 10,
refresh = 1000,
seed = 123,
control = list(adapt_delta = 0.99, max_treedepth = 20),
backend = "rstan"
)Compiling Stan program...
Start sampling
mitchell.brm4 |>
conditional_effects("Day") |>
plot(points = TRUE) |>
_[[1]] +
scale_y_continuous(trans = scales::pseudo_log_trans())mitchell.brm4 |>
conditional_effects("TOD") |>
plot(points = TRUE) |>
_[[1]] +
scale_y_continuous(trans = scales::pseudo_log_trans())9.0.1 stan plots
The brms package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.
Of these, we will focus on:
- stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).
The chains appear well mixed and very similar
- stan_acf (auto-correlation function): plots the auto-correlation between successive MCMC sample lags for each parameter and each chain
There is no evidence of auto-correlation in the MCMC samples
- stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.
All Rhat values are below 1.05, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
Ratios all very high.
9.0.2 pp check
9.0.3 DHARMa residuals
DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components yourself, we can still obtain the simulated residuals from the fitted stan model.
We need to supply:
- simulated (predicted) responses associated with each observation.
- observed values
- fitted (predicted) responses (averaged) associated with each observation
## mitchell.resids <- make_brms_dharma_res(mitchell.brm4, integerResponse = FALSE)
## wrap_elements(~testUniformity(mitchell.resids)) +
## wrap_elements(~plotResiduals(mitchell.resids, form = factor(rep(1, nrow(mitchell_sub))))) +
## wrap_elements(~plotResiduals(mitchell.resids, quantreg = TRUE)) +
## wrap_elements(~testDispersion(mitchell.resids))
## mitchell.resids |> testTemporalAutocorrelation(time = mitchell_sub$Day)
# mitchell.resids |> testTemporalAutocorrelation(time=unique(mitchell_sub$Day))
autocor_check(mitchell_sub, mitchell.brm3, variable = "Day", n.sim = 250)variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
Using as id variables
`geom_smooth()` using formula = 'y ~ x'
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
Using as id variables
`geom_smooth()` using formula = 'y ~ x'
Conclusions:
- the simulated residuals do suggest that there might be a dispersion issue
- it might be worth exploring either zero-inflation, a negative binomial model, or including a observation-level random effect.
mitchell.brm3 |>
augment() |>
## group_by(ID) |>
reframe(ACF = as.numeric(acf(.resid, plot = FALSE)$acf)) |>
## group_by(ID) |>
mutate(N = 1:n()) |>
## slice(1:12) |>
ungroup() |>
group_by(N) |>
summarise(ACF = mean(ACF)) |>
ggplot(aes(y = ACF, x = N)) +
geom_hline(yintercept = 0) +
geom_segment(aes(yend = 0, xend = N))mitchell.brm4 |>
augment() |>
## group_by(ID) |>
reframe(ACF = as.numeric(acf(.resid, plot = FALSE)$acf)) |>
## group_by(ID) |>
mutate(N = 1:n()) |>
## slice(1:12) |>
ungroup() |>
group_by(N) |>
summarise(ACF = mean(ACF)) |>
ggplot(aes(y = ACF, x = N)) +
geom_hline(yintercept = 0) +
geom_segment(aes(yend = 0, xend = N))